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ABSTRACT 

High  order  accurate  centered  flux  approximations  used  in  the  computation  of  numerical 
solutions  to  nonlinear  partial  differential  equations  produce  large  oscillations  in  regions  of 
sharp  transitions.  In  this  paper,  we  present  a  new  class  of  filtering  methods  denoted  by  ENO- 
LS  (Essentially  Nonoscillatory  Least  Squares)  which  constructs  an  upgraded  filtered  solution 
that  is  close  to  the  physically  correct  weak  solution  of  the  original  evolution  equation.  Our 
method  relies  on  the  evaluation  of  a  least  squares  polynomial  approximation  to  oscillatory 
data  using  a  set  of  points  which  is  determined  via  the  ENO  framework. 

Numerical  results  are  given  in  one  and  two  space  dimensions  for  both  scalar  and  systems 
of  hyperbolic  conservation  laws.  Computational  running  time,  efficiency  and  robustness  of 
the  method  are  illustrated  in  various  examples  such  as  Riemann  initial  data  for  both  Burgers’ 
and  Euler’s  equations  of  gas  dynamics.  In  all  standard  cases  the  filtered  solution  appears  to 
converge  numerically  to  the  correct  solution  of  the  original  problem.  Some  interesting  results 
based  on  nonstandard  central  difference  schemes,  which  exactly  preserve  entropy,  and  have 
been  recently  shown  generally  not  to  be  weakly  convergent  to  a  solution  of  the  conservation 
law,  are  also  obtained  using  our  filters. 
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1  Introduction 


Numerical  improvements  in  the  computation  of  high  order  accurate  numerical  solutions  to 
nonlinear  hyperbolic  conservation  laws  have  been  recently  obtained.  Hence,  following  Total 
Variation  Diminishing  (TVD)  schemes,  the  Essentially  Nonoscillatory  (ENO)  method  has 
been  introduced  and  proved  to  be  very  efficient  in  the  computation  of  high  accurate  numerical 
solutions  for  several  types  of  physical  problems  including  Computational  Fluid  Dynamics 
(CFD)  problems  or  front  propagation  using  the  Hamilton-Jacobi  framework.  However,  these 
highfliccurate  methods  use  a  lot  of  computational  time.  For  that  reason,  filtering  methods 

A 

were  developed,  beginning  in  the  late  eighties.  The  first  one,  described  by  B.  Engquist 
and  B.  Sjogreen  in  [1],  uses  simple  TVD  and  conservation  properties  to  correct  nonphysical 
spurious  oscillations  from  one  time  step  to  another.  The  correction  step  consists  in  pushing 
numerical  data  points  up  or  down  to  an  acceptable  level  while  preserving  conservation.  In  [5], 
we  presented  a  new  class  of  filtering  methods  of  any  order  of  accuracy.  Our  method  relies 
on  switching  fluxes  at  locations  in  which  spurious  oscillations  are  detected.  This  method 
was  observed  to  be  very  efficient  and  its  cost  relatively  low  since  high  order  ENO  fluxes  are 
evaluated  at  a  few  points  only-  a  central  difference  method  is  used  most  often. 

In  this  paper,  we  investigate  some  interesting  computational  properties  of  centered  schemes 
after  numerical  oscillations  have  developed  and  propagated  for  some  time.  We  define  a  new 
iass  of  filtering  methods  that  can  be  applied  to  highly  oscillatory  numerical  solutions.  This 
relies  on  the  construction  of  an  ENO  stencil  of  points  ([2,  7,  8])  which  is  fitted  with  high  de¬ 
gree  polynomials  from  a  least  squares  procedure.  Our  numerical  filter  is  capable  of  smoothing 
oscillations  having  large  amplitude  and  high  frequency,  but  without  removing  sharp  singu¬ 
larities  which  are  crucial  components  of  these  solutions.  Furthermore,  the  filtered  solution 
seems  to  retain  the  oscillatory  solution  properties  of  the  unfiltered  schemes  in  some  special 
“entropy  preserving”  cases  as  defined  in  [3]  by  J.  Goodman  and  P.  D.  Lax,  and  in  [6]  by 
J.  Liu  and  D.  Levermore.  We  investigate  some  numerical  examples  using  several  centered 
approximations  in  order  to  illustrate  the  former  property.  The  main  conclusion  indicates 
that  for  standard  central  differences,  our  filtered  solution  always  converges  accurately  to  the 
strong  limit,  whereas  the  predicted  oscillatory  behavior  is  retained  even  after  using  our  filter 
in  the  examples  of  [3,  6]. 

Our  main  test  problems  will  be  inviscid  Burgers’  equation  and  the  inviscid  Euler  equations 
of  gas  dynamics.  We  first  consider  Burgers’  equation: 

«  +  (t).  =  0’  (1) 

with  smooth  initial  condition  U(x,  0)  =  U0(x),  U0  €  Cco(0, 1),  and  periodic  boundary  con¬ 
ditions.  We  will  discuss  the  main  properties  of  the  numerical  solution  obtained  from  some 
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schemes  based  on  the  approximate  fluxes: 


F»i  =  \(U]»  +  Of),  (2) 

FiH  =  &Oj+l  +  Uf))  -  ±(UJM  +  Uf_,),  (3) 

F1+i=  &(0f»  +  m-  UJ-J  +  UUU  +  UU)’  W 

FiH  =  UUJ„  +  UjUj+i  +  Of),  (5) 

Fi+ 1  =  WsOhi-  (6) 


The  fluxes  (2,3,4)  are  just  standard  central  differencing  of  second,  fourth  and  sixth  order 
of  accuracy,  respectively;  while  (5,6)  are  the  interesting  examples  analyzed  in  [6]  and  in  [3], 
respectively.  The  oscillatory  solution  obtained  from  any  of  these  fluxes  is  then  corrected  by 
the  ENO-LS  method,  preserving  the  formal  order  of  accuracy. 

The  Euler  equations  of  gas  dynamics: 

U«  +  (F(U)),  =  0,  (7) 

U(x,0)  =  U0(:r), 

are  to  be  solved  for  t  >  0  and  x  in  some  interval  fl  with  appropriate  boundary  conditions, 
where 

F(U)  =  uU  +  (0  ,p,vp)T 

and  U  =  ( p,q,p ),  p  is  density,  q  is  momentum,  v  is  velocity,  and  p  is  the  pressure.  In  this 
work,  we  use  conventional  second,  fourth  and  sixth  order  central  differencing  with  ENO-LS 
post  processing  applied  to  Euler’s  equations.  See  [4]  for  an  analysis  of  the  oscillatory  Von 
Neumann-Richtmyer  scheme  approximating  Euler’s  equation,  [9]. 

2  The  ENO-LS  Filter,  Algorithms 

The  ENO-LS  method  mimics  the  construction  of  ENO  polynomials  but  without  involving 
the  evolution  equations.  In  short,  we  follow  the  algorithm  just  below: 

Algorithm  2.1  •  1.)  Compute  N  times  the  numerical  solution  of 

Ut  +  A{U )x  =  0 

U{x,  0)=  U0{x) 

i.e  we  let  Vf*  =  Uo(Xj),  and  compute  =  I(V°,NAt),  for  all  j  =  ivhere  I 

is  the  solution  semigroup  operator  that  transforms  the  initial  pointwise  data  Vj°  to 
after  N  iteration  time  steps. 
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•  2.)  Filter  the  numerical  solution  computed  in  step  1.)  by  an  iteration  procedure  similar 
to  Jacobi  or  Gauss-Seidel  elliptic  solvers;  first  let: 

=  V}N 


for  all  j  =  0 then  make  use  of  primitive  variables: 

u?+i  =  'twr*xi, 

2  <=0 

and  finally  construct  a  sequence  {U™  i  }m=o,...,M  s0  that 

J-r  J 

=  E{Um,Um+1), 

JT  2 

where  M  is  defined  from  the  stopping  criterion: 

||wm+1  -  wm||  <£, 


(8) 


(9) 


for  m  —  0, M  —  l,  where  e  is  a  small  parameter  of  order  Axa  and 

w r  =  (io) 

A  Xj 

Finally,  let  VjN  =  W™ ,  for  j  = 

#  3.)  Go  to  step  1.)  unless  t  =  tmax. 

We  notice  that  relation  (10)  and  the  use  of  primitive  variables  (8)  implies  the  conservation 
property  of  the  sequence  {Wj}j= i.e  the  resulting  finite  difference  scheme  is  always  in 
conservation  form.  Moreover,  the  number  of  correction  steps  M  can  be  initially  fixed  so  that 
the  ratio  A.  is  as  large  as  desired.  The  operator  E  is  a  non  trivial  linear  combination  of 
IFF*,},  for  some  j,  as  in  point  Jacobi  or  Gauss  Seidel  method.  In  the  Jacobi  procedure 
we  have: 


and  either 


=  e{u  r+1 

=  E(U?,. 


i )  if  j  varies  from  0  to  ?r,  or 

n  2 

)  in  the  reverse  direction; 


for  the  Gauss-Seidel  method. 
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An  important  property  of  the  ENO-LS  algorithm  comes  from  the  fact  that  the  corrected 
solution  VjN  satisfies  a  conservation  equation.  To  see  that,  we  first  use  the  relation  (10) 
which  can  be  rewriten  as: 


Wm+ 1 

j 


w?+ 


(£'+i(Vj+j-,+,...,£'j+i+.J  -  DJU) -  (Ei-HUi-i- . UH+.J-  U-.) 


and  then  discuss  the  construction  of  the  least  squares  process.  In  (11),  the  pair  of  integers 
(r±,  -Si)  limit  the  width  of  the  stencil  used  in  the  evaluation  of  the  least  squares  polynomial. 
The  appropriate  stencil  is  defined  as  in  the  ENO  algorithms  (refer  to  [7],  [8]).  We  briefly 
indicate  the  main  steps  leading  to  such  a  polynomial:  We  first  compute  the  divided  differences 
table  of  Wm  and  define  the  ENO  stencil  of  points  in  the  region  which  is  the  smoothest  for 
successive  space  derivatives  of  Wm.  We  denote  the  ENO  stencil  in  the  set  (x_,_r, ...,  xJ+3), 
where  r  +  s  —  p  —  1,  and  p  is  the  number  of  data  points  that  we  want  to  take  into  account 
in  the  evaluation  of  the  least  squares  polynomial.  Hence  if  we  denote  this  polynomial  by 

pM(x)  =  £w,(x),  (12) 

i=0 

where  are  the  basis  functions  of  some  polynomial  space  of  degree  q ,  then  the 

unknown  coefficients  Y  =  (Ya,...,Yq)  are  solutions  of  the  linear  system: 

CJ+*  Y  =  Fj+2, 

where  is  a  (q  +  1)  x  (q  +  1)  square  matrix,  F-*+^  is  a  q  +  1  column  vector,  and  both 

can  be  computed  from  the  basis  functions: 

Fk  2=  Hi=-T  Uj+^+iVkiXj+l+i)- 

The  updated  value  follows  from  letting  x  =  xJ+i  on  the  previously  constructed  LS 

polynomial: 

V£V  =  P+Htni)  = 


The  global  conservation  feature  follows  provided  that  we  write: 

S  1 

TTm+ 1  -  V'  ry-7+2  ffm+l 

uj+k  ~  ^ 

2  t=-r  2 

where  /  =  0  or  1  depending  on  whether  the  Jacobi  or  the  Gauss-Seidel  method  is  used,  and 


=  £  £  D\+k 


j+k+t> 


t=0  k= 0 


t 


where  DJ+  2  =  ( C J+2)  E 

j+i 

Note  first  that  the  coefficients  {a  ?  }t=_r . 3  form  a  sequence  of  bounded  real  numbers, 

and  second,  that  the  basis  functions  (<^>o>  •  ■•Tq)  can  be  appropriately  chosen  as  a  sequence  of 
Chebyshev  polynomials  or  some  other  set  of  orthogonal  polynomials.  The  last  choice  permits 
us  to  compute  directly  the  coefficients  without  inverting  the  matrix  CJ+ 5: 


a 


J  +  2 


j+h+t 


i=o 
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This  yields  a  fast  algorithm  since  matrix  inversions  are  no  longer  needed. 

We  now  present  an  extension  of  algorithm  2.1  to  two  space  dimensions.  The  simplest 
possible  extension  would  be  to  apply  the  previous  algorithm  in  two  sweeps.  The  first  one 
will  freeze  one  coordinate  and  correct  the  oscillatory  solution  with  respect  to  the  other  free 
variable.  The  second  one  will  simply  reverse  the  role  of  each  variable.  This  method,  while 
simple,  has  difficulties  near  curved  shocks.  We  shall  use  instead  a  fully  two  dimensional  ENO- 
LS  filtering  algorithm.  The  latter  will  provide  the  construction  of  least  squares  polynomials 
in  two  space  dimensions  using  a  set  of  points  which  is  chosen  as  the  intersection  of  one 
dimensional  ENO  stencils  of  data  points  in  each  separate  direction.  The  algorithm  below 
describes  our  procedure. 


Algorithm  2.2  •  1.)  Compute  N  times  the  numerical  solution  of 

Ut  +  A(U)X  +  B(U)y  =  0 

U(x,y,  0)  =  U0{x,y) 

using  a  very  simple  numerical  method  and  then  let  =  Uo(xj,y,),  and  V/j  = 
I(V°,NAt),  for  all  j  =  1, ..., Ti\,  and  i  =  l,...,n2,  where  I  is  the  solution  semigroup 
operator  that  we  have  already  encountered  in  previous  algorithm. 

•  2.)  Let  Wfj  =  Vfj ,  and  filter  M  times  the  primitive  variable 

1,)+1  =  E  W^Ax^Aj,,,  (13) 

l,k=0 

by  a  Jacobi  or  point  Gauss-Seidel  iteration  procedure,  i.e  perform: 

u”y\\+>  =  E(um,um+l), 

for  positive  integers  m.  Iterate  as  long  as  ||Wm+1  —  Wm\\  <  e,  in  =  0,  —  1,  and 

finally  let  the  filtered  solution:  V/V  =  ,  for  j  —  l,...,«i,  1  —  l,...,ri2. 

The  construction  of  the  two  dimensional  least  squares  polynomial  is  as  follows: 
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•  2.1)  At  the  location  (xi+i,yj+i),  compute  the  ENO  stencils  of  points 
{(•£«— 2/j ) » •••»  (■^t+i+sij and  {(:£»»  Ty ))•••?  (*£ii  J/i+i+sy)}>  /or  vx  4-  sx  4*  1  —  ri\ 
and  ry  +  sy  +  1  =  n2,  where  nj  and  n2  are  the  preset  maximum  number  of  points  along 
the  x  and  y  axis.  These  sets  of  points  form  two  segments  crossing  at  (x{,yj).  Then , 
define  the  x  and  y  ENO  stencils  of  points  along  these  y  and  x  segments,  respectively. 
The  two  dimensional  ENO  stencil  is  taken  as  the  intersection  of  the  union  of  the  pre¬ 
defined  x  and  y  one  dimensional  ENO  stencils  (refer  to  figure  (1)).  The  least  squares 
polynomial  is  then  simply  defined  on  this  two  dimensional  ENO  stencil  and  is  set  to 
p('+7’i+5)(x,y)  =  Y)l=oYkTk(x,y)-  Again,  the  unknowns  coefficients  Yk,  k  =  1  ,...,p 
are  computed  by  solving  the  linear  system  Cb+jj+lly  =  Fb+id+j/ 


.  2.11)  Let  V~ |‘J+,  = 

•  2.3)  Recover  the  conservative  solution 

uR'i  +  U?t\  .  -  uz tv  ,  -  ur t*. 


WmAl  = 


•+5J+5 


_* _ 2 


'+  2  J  2 


i-kd+i 


Ax  Ay 


for  m  =  0, M  —  1. 


(14) 


•  3.)  Go  to  step  1.)  unless  t  =  tmax. 


An  example  of  two  dimensional  ENO  stencil  is  given  in  figure  1.  The  main  interesting 
feature  of  such  construction  is  based  on  the  localization  of  the  least  oscillatory  part  of  the 
solution  within  the  much  larger  rectangle  Xj+Sx]  x  [yi_r„>y«+*y]- 

In  section  4,  we  investigate  the  two  dimensional  Burgers’  equation  and  study  numerical 
propagation  of  a  shock  along  the  radial  axis.  With  centered  fluxes,  some  spurious  numerical 
oscillations  propagate  in  the  direction  of  the  flow;  however  the  two  dimensional  ENO-LS 
filtering  method  is  able  to  filter  all  these  oscillations  while  still  giving  the  correct  location  of 
the  curved  shock  wave. 

Next,  we  consider  hyperbolic  systems  of  conservation  laws  exemplified  by  the  inviscid 
Eulers’  equations  of  compressible  gas  dynamics: 

Ut  +  F(U)X  =  0, 

U(x,0)  —  U0{x) 

for  which  there  exists  a  complete  set  of  real  eigenvectors  and  eigenvalues;  i.e  NF(U)  = 
P~lAP,  where  A  is  a  diagonal  matrix  with  entries  Aj  <  A2...  <  A j,  and  P  is  a  matrix  whose 
columns  define  a  complete  set  of  eigenvectors  of  the  system.  Note  that  the  eigenvalues  can  be 
multiple,  which  is  the  case  for  the  Euler  equations  of  gas  dynamics  in  two  space  dimensions. 
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Figure  1:  2D  ENO-LS  stencil. 


Using  the  eigenvector  decomposition,  a  field  by  field  approach  is  used,  i.e  the  fluxes  will  be 
corrected  in  each  fields  (our  filter  sometimes  failed  to  remove  oscillations  near  strong  shocks 
when  applied  to  the  conserved  quantities).  Briefly,  we  proceed  as  follows: 


Algorithm  2.3  (Field  by  field  ENO-LS  method)  •  2-1.)  Let  Wm  =  {Wf, ...,  W™)T , 
and  define  the  primitive  variable  V™  i  =  £<=i  W™Az. 

•  2-2.)  Compute  the  left  and  right  eigenvectors  L"V  ^ for  k  =  l,...,d  (d  is  the 
number  of  equations)  of  the  Roe  matrix  (see  e.g.  [2])  and  decompose 

the  primitive  vector  V™(_1  along  each  characteristic  field: 


m,(k)  _ 

J+l 


^  '  2 
m»(*) 


Ed 

fe=l  <7,+  i 


,R’ 


».(*) 


/or  which  we  have  frozen  the  index  j  in  order  to  get  the  same  decomposition  for  all 
neighboring  points  for  l(k)  =  —  r(k), ...,  +s(fc)  involved  in  the  calculation  of 

the  least  squares  polynomial  which  is  constructed  in  step  2-3.). 

•  2-3.)  Select  an  ENO  stencil  of  p  points,  i.e  {xj-T(k), ...,  £_,+,(*)},  for  r(k)  +  s(k)  +  1  —  p; 
and  then  define  the  least  squares  polynomial  of  degree  q  so  that: 


j+k+,(k)' 
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•  2-4-)  Transform  back  the  filtered  field  by  field  solution  to  the  primitive  vector: 


v74'  =  £  CiMR7;f’ 


k=l 


and  finally  recover  the  desired  physical  variables: 

■yrn+l  _  -y-m+l 


m+1  _  i 


w™+1  = 


i- 


Ax 


(15) 


•  2-5.)  Iterate  until  the  stopping  criterion 


||Wm+l  -  Wm||  <  £, 

for  m  =  0, M  —  1,  is  reached;  and  finally  let  =  Wj^. 

Note  that  this  algorithm  can  be  extended  to  two  space  dimensions  by  correcting  separately 
oscillatory  fields  involved  in  the  x  and  y  fluxes,  respectively.  Moreover,  this  algorithm  does 
not  make  use  of  the  evolution  equations  but  does  use  the  eigenfunction  decomposition  in 
order  to  track  efficiently  the  propagation  of  spurious  oscillations. 

To  conclude  this  section,  we  indicate  that  we  mainly  supposed  that  the  numerical  oscil¬ 
lations  always  propagate  with  the  flow  speed  along  local  characteristic  fields  and  that  the 
amplitudes  of  such  oscillations  are  not  too  large  so  that  the  oscillatory  solution  does  not  be¬ 
come  unbounded,  e.g.,  negative  density  and/or  pressure  is  not  allowed.  In  all  our  numerical 
experiments,  we  had  to  turn  the  filter  on  not  only  to  recover  an  acceptable  final  solution  but 
also  to  reduce  the  amplitude  and  frequency  of  spurious  oscillations  during  the  calculations. 
Hence,  we  usually  preset  the  value  of  the  ratio  ^  in  the  numerical  experiments  just  after 
singularities  have  developed. 


3  Numerical  Convergence  Study 

In  this  section,  we  investigate  the  numerical  convergence  of  the  approximate  solution  com¬ 
puted  via  the  ENO-LS  algorithm  given  an  initial  oscillatory  solution  which  has  been  eval¬ 
uated  from  one  of  the  standard  centered  fluxes  (2),  (3),  (4);  or  from  one  of  the  ’’entropy 
conserving”  fluxes  (5)  or  (6).  As  test  problem,  we  study  first  the  numerical  evolution  of 
Burgers’  equation  (1)  in  one  space  dimension: 

U2 

tt  +  (y)*  =  o, 

with  initial  condition  U(x,  0)  =  sin  27rx,  in  the  domain  [0,1],  and  extend  the  solution  by 
periodicity  outside  0  and  I.  A  shock  wave  develops  at  t  =  ~  at  the  point  x  = 
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We  first  display  in  figure  (2)  the  solution  using  100  cells  and  a  CFL  condition  of  | 
for  a  second  order  centered  difference  (CD)  scheme  using  the  flux  (2),  and  then  correct 
that  solution  by  applying  7  iterations  of  the  ENO-LS  algorithm.  Note  first  that  the  exact 
solution  is  perfectly  recovered  with  exact  location  of  the  singularity  because  of  the  use 
of  primitive  variables  in  algorithm  2.1;  second,  these  results  are  obtained  for  the  set  of 
coefficients  (p,  q)  =  (13,2),  which  implies  that  the  corrected  solution  is  only  secjnd  order 
accurate  in  smooth  transition  areas;  third,  the  upper  right  plot  displays  the  corrected  solution 
when  both  (3,2)CD  and  (3,2)ENO  schemes  are  sequentially  used.  In  that  case,  we  correct  by 
5  ENO  iterations  after  every  25  CD  steps.  This  method  for  filtering  oscillations  is  also  very 
efficient  for  one  dimensional  problems  (refer  to  table  1);  however,  many  more  ENO  iterations 
are  needed  in  two  dimensions.  Basically,  10  ENO  iterations  are  needed  after  every  10  CD 
steps  in  order  to  recover  an  accurate  solution  for  the  two  dimensional  Burgers’  example  4.3, 
which  is  a  quite  expensive  technique  compared  to  the  overall  cost  induced  by  the  ENO-LS 
method. 

We  shall  discuss  the  results  obtained  with  larger  values  of  q  in  section  4  in  which  we 
investigate  numerical  order  of  accuracy  of  the  filtered  solutions  from  the  ENO-LS  algorithm. 
Also,  a  similar  study  is  investigated  as  the  number  of  evaluation  points  p  increases. 

In  figure  (3),  we  plot  the  numerical  solution  before  and  after  the  filter  steps  when 
n  =  1000.  Again,  convergence  to  the  physical  solution  is  reached  within  10  Gauss  Sei¬ 
del  iterations.  Note  that  this  number  increases  by  a  factor  of  nearly  2  when  Jacobi  iterations 
are  used. 

In  figures  (4)  and  (5),  we  show  the  filtered  solutions  which  were  initially  computed  using 
the  standard  (3,4)CD  and  (3,6)CD  schemes  (fluxes  (3)  and  4)),  respectively.  Note  that  the 
corrected  solution  is  well  reconstructed  in  smooth  regions  while  a  smearing  of  the  shock  over 
about  10  cells  is  observed.  This  smearing  can  be  primarily  explained  from  the  high  (fourth) 
degree  least  squares  polynomials  ( q  =  4)  used  in  those  experiments. 

Our  second  numerical  test  is  devoted  to  showing  that  the  filtered  solution  computed  via 
the  numerical  flux  (5): 


FJ+l=  +  IW+.  +  Uj). 

does  not  converge  to  the  physical  solution.  It  was  observed  in  [6]  that  this  scheme  preserves 
mass  and  entropy  (which  in  this  case  is  taken  to  be  \U2).  It  was  also  shown  there  that  the 
numerical  solution  does  not  converge  to  a  weak  solution  of  the  original  problem.  Basically, 
as  the  mesh  size  tends  to  zero,  some  spurious  oscillations  are  still  visible  near  the  shock  and 
cannot  be  removed.  We  ran  the  last  two  previous  experiments  but  used  the  approximate  flux 
(5).  We  plot  the  numerical  results  in  the  figures  (6)  and  (7)  when  n  =  100  and  n  =  1000, 
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respectively.  For  100  cells  (figure  (6)),  the  oscillations  near  the  shock  are  all  smoothed  out; 
however,  the  location  of  the  shock  is  smeared  over  about  10  cells.  This  results  were  obtained 
for  the  same  set  of  coefficients  (p,  q)  =  (13,2).  Indeed,  the  solution  is  not  well  reconstructed. 
As  the  number  of  cells  increases,  the  filtered  solution  is  well  reconstructed  except  near  the 
shock.  Numerical  oscillations  are  still  visible  and  cannot  be  removed.  Note  that  the  results 
visualized  in  figure  (7)  are  given  after  16  Gauss  Seidel  iterations  which  is  much  more  than 
we  needed  when  standard  centered  differences  are  taken. 

Our  final  Burgers’  equation  test  deals  with  a  similar  convergence  failure  property  to  the 
correct  physical  solution  when  the  flux  approximation  (6)  is  implemented.  The  numerical 
flux  is: 

This  scheme  again  conserves  both  mass  and  entropy  (this  time  the  entropy  is  taken  to  be 
log Uj),  and  it  was  shown  in  [3]  that  the  numerical  solution  does  not  converge  to  a  weak 
limit  of  the  original  problem  as  the  stepsize  Ax  tends  to  zero.  In  the  numerical  experiments, 
we  noticed  that  the  amplitude  of  the  oscillations  grew  very  fast  and  was  rapidly  becoming 
unbounded.  The  results  are  plotted  in  the  figures  (8)  and  (9)  for  n  —  100  and  n  =  1000, 
respectively.  Note  that,  in  the  case  n  =  100,  the  filtered  solution  is  again  not  very  well 
reconstructed,  and  for  1000  cells,  some  oscillations  are  still  visible  after  20  ENO-LS  iterations 
on  the  right  of  the  shock.  Moreover  these  oscillations  could  not  be  removed,  even  after  many 
additional  filtering  iterations. 

Thus  our  techniques  have  been  observed  to  construct  converging  sequences  of  filtered 
numerical  solutions  towards  the  expected  solution  given  initial  oscillatory  data  that  has 
been  computed  from  a  central  difference  flux  approximation.  Our  hope  is  now  to  show  that 
our  method  is  not  only  robust  but  is  also  fast.  This  is  the  topic  of  last  section. 

4  Time  Efficiency  of  ENO-LS  Algorithms 

In  this  section,  we  want  to  test  the  ENO-LS  filtering  method  for  several  test  problems 
involving  nonlinear  hyperbolic  systems  of  conservation  laws.  We  will  focus  our  attention 
on  comparing  precisely  the  CPU  times  of  the  ENO-LS  method  versus  more  classical  and 
filtered  methods.  Among  them,  we  will  consider  the  straightforward  central  difference  (CD) 
method,  the  expensive  ENO  (Essentially  Nonoscillatory)  technique  [7,  8],  and  our  FM  scheme 
(Filtering  Method)  [5].  We  will  run  three  examples  for  ID  and  2D  Burgers’  equation,  and 
for  Eulers’  equations  of  compressible  gas  dynamics. 
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Type  of  Scheme 

Order  of  Accuracy 

CPU  time  xlO"2 

#  of  corrections,  comments 

CD 

(3,2) 

0.56 

CD  =  Centered  Differences 

CD 

(3,4) 

0.72 

- 

ENO-RF 

(3,2) 

1.11 

RF  =  Roe  Fix,  Entropy  fix 

ENO-RF 

(3,4) 

1.68 

- 

FM 

(3,2) 

0.76 

FM  =  Filtering  Method,  3%  Corrections 

FM 

(3,4) 

1.07 

4%  Corrections 

FM 

(3,2) 

1.40 

100%  Corrections  =  full  ENO 

FM 

(3,4) 

2.04 

100%  Corrections 

CD+ENO  (40,5,20) 

(3,2) 

0.74 

40CD,  5ENO,  20  final  ENO 

CD+ENO  (40,1,10) 

(3,2) 

0.58 

40CD,  1ENO,  10  final  ENO 

CD+ENO  (20,7,10) 

(3,4) 

0.99 

20CD,  7ENO,  10  final  ENO 

ENO-LS  (7,2) 

(2) 

1.84 

(p,  q)  —  (7,2),  LU  Inversion 

ENO-LS  (7,2) 

(2) 

0.65 

Orthogonal  polynomials 

ENO-LS  (7,3) 

(3) 

2.5 

LU  Inversion 

ENO-LS  (7,3) 

(3) 

0.89 

Orthogonal  polynomials 

Table  1:  CPU  times  of  CD,  ENO,  FM  and  ENO-LS  methods. 


4.1  ID  Burgers’ 


We  first  compare  the  time  efficiency  of  several  numerical  schemes  for  the  ID  Burgers’  equation 
(1).  In  table  (1),  we  show  the  average  computational  time  per  iteration  for  100  mesh  points 
for  the  CD,  ENO,  and  FM  schemes  with  several  correction  angles  (see  [5]),  and  for  various 
combinations,  performing  alternatively  some  centered  difference  and  some  ENO  steps.  The 
notation  CD+ENO  (40,5,20)  simply  means  that  the  calculation  starts  with  40  CD  steps, 
followed  by  5  full  ENO  iterations,  and  back  to  the  centered  scheme:  the  last  20  steps  of 
the  calculations  are  finally  performed  by  the  ENO  method  in  order  to  recover  an  acceptable 
solution.  Note  that  all  coefficients  displayed  in  this  table  are  tuned  so  that  the  numerical 
solution  is  high  order  accurate  in  smooth  regions  and  no  spurious  oscillations  are  detected. 

The  contents  of  this  table  need  a  few  comments.  First  of  all,  the  fastest  algorithm  is  the 
one  based  on  central  differences.  This  is  indeed  not  surprising  since  only  one  Fortran  instruc¬ 
tion  is  needed  in  the  coding  of  the  approximate  flux.  Second,  postprocessing  a  numerical 
solution  computed  from  the  CD  scheme  by  an  ENO  method  can  be  very  efficient  for  lower 
orders.  For  a  second  order  method,  only  one  ENO  iteration  after  every  40  CD  steps  has  to  be 
implemented  in  order  to  reduce  sufficiently  the  amplitude  of  oscillations.  Note  however  that 
for  the  fourth  order  method,  7  ENO  iterations  were  needed  every  after  only  20  CD  steps. 
In  fact,  if  these  spurious  oscillations  are  not  regularly  cut  off,  the  final  ENO  iterations  may 
not  recover  the  correct  numerical  solution.  Third,  the  ENO-LS  filtering  method  is  the  most 
costly  when  full  LU  inversion  of  the  CJ+2  matrix  is  performed.  However,  when  orthogonal 
basis  functions  are  introduced,  the  ENO-LS  method  is  competitive  with  respect  to  the  fast 
CD  scheme.  Note  moreover  that  ENO-LS  correction  steps  have  to  be  performed  at  a  few 
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Type  of  Scheme 

(q.p) 

#  of  iterations 

L 1  and  L°°  orders. 

CD  (3,2) 

(2,5) 

12 

1.81  and  1.63 

CD  (3,2) 

(2,9) 

9 

1.75  and  1.50 

CD  (3,2) 

(2,13) 

7 

2.19  and  1.9 

CD  (3,4) 

(3,5) 

15 

1 .59  and  1 .40 

CD  (3,4) 

(3,9) 

12 

2.41  and  2.33 

CD  (3,4) 

(3,13) 

7 

2.80  and  2.79 

Table  2:  Local  Ll  and  Z,°°  order  of  accuracy. 


times  only.  Finally,  after  running  many  experiments,  we  noticed  that  if  the  ratio  ^  becomes 
too  large,  then  shocks  have  a  tendency  to  spread  over  a  large  number  of  cells.  Again,  there 
is  a  compromise  that  needs  to  be  reached  for  fast  convergence;  this  depends  on  the  large 
value  of  the  ratio  and  the  approximation  of  the  shape  near  shocks  for  which  p  needs  to 
be  slightly  smaller.  Note  that  in  most  of  experiments,  the  ratio  ^  €  [4,6]  was  optimal. 

Now,  we  want  to  measure  the  order  of  accuracy  of  the  ENO-LS  solution.  To  do  so,  we 
measure  from  computations  the  L 1  and  L°°  errors  in  the  slabs  [0.10,0.24]  and  [0.26,0.30]. 
Numerical  orders  are  shown  in  table  2. 

Several  comments  about  these  results  are  now  discussed.  First  of  all,  as  the  number  of 
evaluation  points  increases,  the  better  the  quality  of  the  results  and  the  faster  convergence 
is  reached.  Indeed,  we  have  to  pay  the  price  of  higher  computations  which  are  required  to 
construct  the  coefficients  involved  in  the  matrix  and  in  the  vector  FJ+* .  On  the  other 

hand,  faster  calculations  can  be  obtained  provided  that  the  values  of  p  and  q  diifer  only 
slightly,  i.e  p  =  q  +  /,  /  =  1,2,3,...,  but  local  accuracy  becomes  obviously  poor.  Again,  the 
optimal  ratio  seems  to  belong  to  the  interval  [4,6]. 

4.2  Example  2.  Euler  Equations  of  Gas  Dynamics 

The  second  test  problem  is  devoted  to  the  Euler  equations  of  gas  dynamics  in  one  space 
dimension.  We  consider  the  initial  condition  given  in  example  8  of  [8],  that  is: 

p  =  3.857143,  q  =  2.629369,  p  =  10.3333333  when  x  <  -4 
p  =  1  +  £su\5x,q  =  0,p  =  1.  when  x  >  — 1 


When  e  =  0,  a  pure  Mach  3  shock  is  moving  to  the  right  from  the  initial  discontinuity 
x  =  4.  When  £  =  0.2,  we  not  only  have  a  Mach  3  shock  propagating  to  the  right,  but  have  as 
well  a  succession  of  weaker  rarefaction  and  shock  waves  propagating  to  the  left.  Numerical 
results  for  the  ENO  and  FM  methods  of  high  order  of  accuracy  can  be  found  in  [8]  and  in  [5], 
respectively.  We  ran  the  same  problem  using  240  CD  iterations  and  then  plotted  the  results 
in  figure  (10).  Next,  we  correct  this  highly  oscillatory  numerical  solution  by  performing  7 
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Figure  10:  (3,2)CD  scheme. 

ENO-LS  correction  steps.  In  this  experiment,  we  use  (p,  q)  =  (13,3).  The  filtered  solution 
is  visualized  in  figure  (11).  Note  that  the  pressure,  velocity,  and  entropy  are  quite  well 
reconstructed,  whereas  the  density  is  not  perfectly  recovered  near  the  strong  shock  for  which 
some  physical  oscillations  should  remain  (  refer  to  [10,  11]).  However,  we  should  note  the 
remarkable  improvement  obtained  from  the  solutions  displayed  in  figures  (10)  and  (11). 

In  figures  (12)  and  (13),  we  use  in  sequence  50  CD  steps  and  only  one  ENO  iteration 
for  second  and  fourth  order  methods  before  correcting  the  final  oscillatory  solution.  The 
oscillatory  solution  is  now  postprocessed  by  3  ENO-LS  iterations.  The  filtered  numerical 
results  now  look  fairly  similar  to  those  shown  in  [8,  5]. 

4.3  Example  3.  2D  Burgers’  Equation 

The  last  example  is  devoted  to  the  two  dimensional  Burgers’  equation  to  be  solved  in  the 
square  domain  [— 1, 1]  x  [— 1, 1]  with  initial  condition  Uo{x,y)  =  sin27rr,  where  r  =  y/  x2  +  p2, 
for  r  <  and  Uo{x,y)  =  0  outside  the  disc  r  =  In  figure  (14),  we  visualize  the  solutions 
obtained  by  the  (3,2)CD  scheme,  followed  by  4  iterations  of  ENO-LS  correction  steps.  In 
that  experiment,  we  use  (pi,py,<7)  =  (6,6,2),  so  that  the  local  rectangle  in  which  the  two 
dimensional  ENO-LS  stencil  of  points  is  taken  contains  36  mesh  points.  Note  that  a  speed  up 
factor  of  nearly  1.5  is  obtained  by  using  the  two  dimensional  ENO-LS  method  instead  of  the 
one  dimensional  splitting  version.  In  this  numerical  example,  the  amplitude  of  the  numerical 
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Figure  11:  (3,2)CD  scheme  +  7  ENO-LS  method  (p,q)=(13,3). 

oscillations  which  propagate  radially  away  from  the  initial  shock  location  was  approximately 
|  of  the  strength  of  the  shock,  yet  the  ENO-LS  method  recovered  the  nonoscillatory  accurate 
solution  quite  well. 


5  Concluding  Remarks 

The  notion  of  using  an  essentially  nonoscillatory  (ENO)  least  squares  filter  to  remove  spu¬ 
rious  oscillations  of  data  generated  by  numerical  overshoot  appears  promising.  Tuning  of 
parameters  is,  however,  still  required.  Future  work  will  deal  with  this  issue. 
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-PRESSURE- 


(3 ,2)CD+ENO  (50,1)  +  5  final  (ENO-LS)  (p,q)=(13,3)  - 


Figure  12:  (3,2)(CD,ENO)  schemes  +  3  ENO-LS  method. 


(3 ,4)CD+ENO  (50,1)  +  3  final  (ENO-LS)  (p,q)=(13,4)  - 


Figure  13:  (3,4)(CD,ENO)  schemes  -f  3  ENO-LS  method. 


-  3  2D  (ENO-LS)  Jacobi  steps  - 

Figure  14:  (3,2)CD  scheme  +  2D  ENO-LS  method. 
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